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The energy injection and losses in the Monte Carlo simulations 
of a diffusive shock 
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Abstract. Although diffusive shock acceleration (DSA) could be simulated by some well- 
established models, the assumption of the injection rate from the thermal particles to 
the superthermal population is still a contentious problem. But in the self-consistent Monte 
Carlo simulations, because of the prescribed scattering law instead of the assumption of 
the injected function, hence particle injection rate is intrinsically defined by the prescribed 
scattering law. We expect to examine the correlation of the energy injection with the 
prescribed multiple scattering angular distributions. According to the Rankine-Hugoniot 
conditions, the energy injection and the losses in the simulation system can directly de- 
cide the shock energy spectrum slope. By the simulations performed with multiple scat- 
tering law in the dynamical Monte Carlo model, the energy injection and energy loss func- 
tions are obtained. As results, the case applying anisotropic scattering law produce a small 
energy injection and large energy losses leading to a soft shock energy spectrum, the case 
applying isotropic scattering law produce a large energy injection and small energy losses 
leading to a hard shock energy spectrum. 



1. Introduction 

The gradual solar energetic particles with a power-law 
energy spectrum are generally thought to be accelerated by 
the first-Fermi acceleration mechanism at the interplanetary 
shocks (IPs) [Axford et ai, 1977; Krymsky, 1977; Bell, 1978; 
Blandford and Ostriker, 1978]. It is well known that the dif- 
fusive shock accelerated the particles efficiently by the accel- 
erated particles scattering off the instability of Alfven waves 
which are generated by the accelerated particles themselves 
[Lagage and Cesarsky, 1983; Gosling et. al. , 1981; Cane 
et.al. , 1990; Lee and Ryan, 1986; Pelletier et.al., 2006; Li 
et at, 2009]. The diffusive shock acceleration (DSA) is so 
efficient that the back-reaction of the accelerated particles 
on the shock dynamics cannot be neglected. So the theo- 
retical challenge is how to efficiently model the full shock 
dynamics [Capnoli et. al, 2010; Zank, 2000; Li et al, 
2003; Lee, 2005]. To efficiently model the shock dynam- 
ics and the particles' acceleration, there are largely three 
basic approaches: stationary Monte Carlo simulations, fully 
numerical simulations, and semi-analytic solutions. In the 
stationary Monte Carlo simulations, the particle popula- 
tion with a prescribed scattering law is calculated based on 
the particle-in-cell (PIC) techniques [Ellison et al., 1996; 
Vladimirov et al., 2006]. In the fully numerical simulations, 
a time-dependent diffusion-convection equation for the CR 
transport is solved with coupled gas dynamics conservation 
laws [Kang and Jones, 2007; Zirakashvili and Aharonian, 
2010]. In the semi-analytic approach, the stationary or 
quasi-stationary diffusion-convection equations coupled to 
the gas dynamical equations are solved [Blast et. al., 2007; 
Malkov et. al, 2000]. Since the velocity distribution of su- 
perthermal particles in the Maxwellian tail is not isotropic 
in the shock frame, the diffusion-convection equation can- 
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not directly follow the injection from the non-diffusive ther- 
mal pool into the diffusive CR population. So consider- 
ing both the quasi-stationary analytic models and the time- 
dependent numerical models, the injection of particles into 
the acceleration mechanism is based on an assumption of 
the transparency function for thermal leakage [Blasi et. al., 
2005; Kang and Jones, 2007; Vainio and Laitinen, 2007] in 
priori. Thus, the dynamical Monte Carlo simulations based 
on the PIC techniques are expected to model the shock 
dynamics time-dependently and also can eliminate the sus- 
picion arising from the assumption of the injection [Knerr 
et. al., 1996; Wang and Yan, 2011]. In plasma simulation 
(Monte Carlo model and hybrid model), since the proton' 
mass is very larger than the electron' mass, the total plasma 
can be treated as one species of proton fluid with a massless 
electronic fluid which just balance the electric charge state 
for maintaining a neutral fluid [Leroy et al., 1982]. There 
is no distinction between thermal and non-thermal parti- 
cles, hence particle injection is intrinsically defined by the 
prescribed scattering properties, and so it is not controlled 
with a free parameter [Caprioli et. al., 2010]. 

Actually, Wang and Yan [2011] have extended the dy- 
namical Monte Carlo models invoking multiple scattering 
angular distributions. Unlike the previous KJE[_R'nerr et. 
al., 1996] dynamical Monte Carlo models invoking a purely 
isotropic scattering angular distribution, this multiple scat- 
tering law allow the particles are scattered by angles dis- 
tributed with Gaussian functions. According to the simu- 
lations using the extended multiple scattering angular dis- 
tributions, a series of similar energy spectrums with a little 
difference with respect of the power-law tail are obtained. 
And the results show that the energy spectral index is ef- 
fected by the prescribed scattering law. Specifically, the to- 
tal shock's energy spectral index is less than one and shows 
an increasing function of the dispersion of the scattering an- 
gular distribution, but the subshock's energy spectral index 
is more than one and shows a decreasing function of the 
dispersion of the scattering angular distribution. 

In an effort to research why the multiple scattering an- 
gular distributions can produce the difference of the energy 
spectral index, it is necessary to analyze the energy injec- 
tion and the energy losses in the entire simulation system. 
Because the energy injection and losses are important fac- 
tors for deciding the acceleration efficiency and the energy 
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spectrum slope owing to the Rankine-Hugoniot relationship 
based on the energy conservational law. 

However, in the Monte Carlo simulation, the particle in- 
jection and the energy loss processes arc treated in natural, 
self-consistent manner and decided by the prescribed scat- 
tering law. In order to obtain the complete energy injection 
and loss real-time functions in the entire simulation system, 
we perform the simulations by the multiple scattering law 
considering an improved simulation system. In this new sim- 
ulation system, a radial reflective boundary(RRB) is set for 
preventing the energy losses via the radial diffusion. Un- 
der these scenarios, the performed simulation cases consist 
of four specific standard deviation values of the Gaussian 
distribution function. 

In Section 2, the basic simulation method is introduced 
with respect to the Gaussian scattering angular distributions 
for obtaining the energy injection and loss functions of time 
in each case. In Section 3, we present the energy analysis 
for all cases with four types of scattering angle distributions. 
Section 4 includes a summary and the conclusions. 

2. Method 

The Monte Carlo model is a general model, although it is 
considerably expensive computationally, and it is important 
in many applications to include the dynamical effects of non- 
linear DSA in simulations. Since the prescribed scattering 
law can replace the electromagnetical field calculation which 
is used in hybrid simulations [Giacalone, 2004; Winske and 
Omidi, 2011], we assume that the individual particle scat- 
ters elastically off the background scattering centers with 
the scattering angles according to a Gaussian distribution 
in the local frame. And the particle's mean free path is 
proportional to the local velocities in its local frame with 

A = Vl • r. (1) 

Where, t is the average scattering time. Under the pre- 
scribed scattering law, the injection is purely correlated with 
those particles from the "thermal pool" in the downstream 
region become into the superthermal particles [Ellison et al., 
2005]. 

In these simulations, the entire shock is simulated in one- 
dimensional box as shown in Figure 1, the initial continually 
inflow enter into the box from the left boundary with a su- 
personic bulk velocity (Uo), a stationary reflective wall at 
the right boundary of the box act to form a piston shock 
moving from right to left. After a certain time, a steady 
compression region (i.e. downstream region) will be formed 
in front of the reflective wall. The bulk velocity in down- 
stream region is become to zero, since the particles dissipate 
in the downstream region and their large translational en- 
ergy is converted into isotropic, random energy. To model 
the flnite size of system and the lack of sufficient scatter- 
ing far upstream to turn particles around [Mitchell et. al., 
1983], the presented simulation includes the escape of the 
energetic particles at an upstream "free escape boundary" 
(FEB). This FEB moves with the shock front at a shock ve- 
locity {Vsh) and remains a constant distance in front of the 
shock position (i.e. Xfeb~90). This distance is enough 
large for majority of the injected particles diffuse between 
the forcshock region and the downstream region. The size 
of the forcshock is the distance from the shock to the FEB 
and thus sets a limit on the maximum energy a particle 
can obtain. Since the injected particles cross the shock and 
diffuse upstream, they negatively contribute to the bulk ve- 
locity, and the bulk velocity become smaller and smaller 
from the FEB to the shock position. Holding the length of 
the forcshock region constant eventually (when enough time 



has elapsed to create a larger rmmber of accelerated par- 
ticles) produces a steady state with respect to the amount 
of the energy entering and exiting the system from the up- 
stream region. In addition, we set the radial boundary as 
Ry^z = 50 for preventing particle's perpendicular diffusion 
to the infinity. Simultaneously, The radial reflective bound- 
ary can ensure the particles have efficient diffusive processes 
in the one-dimensional system along x direction. So we are 
able to compare the difference of the energy injection and 
losses obtained from each case. We can further investigate 
the possibility that the cases applying anisotropic scattering 
angular distribution would produce a different acceleration 
efficiency compared with the case applying an isotropic scat- 
tering angular distributions. So an anisotropic scattering 
law in the theory of the CR-diffusion is also needed [Bell, 
2004]. 

According to particle-in-cell (PIC) techniques [Forslund, 
1985; Spitkovsky, 2003; Nishikawa et. al, 2008], the total 
box length in this simulation system is Xmax~300, and it 
is divided up into na;=600 grids. The initial number of 
particles in each grid is no =650. In addition, we use a 
flux-weighted inflow to ensure the particles entering into 
the box with the same density flow in upstream with the 
time. This inflow in "preinflow box" (PIB) is put in the left 
boundary of the simulation box. The total simulation time 
Tmax = 2400, and it is divided into the number of time steps 
Nt = 72000 with a time step dt = 1/30. The size of the 
FEB distant from the shock front is set as Xfet = 90. The 
radii of the radial reflective boundary is set as Ry,z = 50. 
These simulation codes consist of the three substeps. (i) In- 
dividual particles move along the x, y, and z ikxis with their 
local velocities in each component, respectively. 

x = xo+Vx -t y = yo+Vy -t z = zo + Vz ■ t (2) 

Since the magnetic field Bq is parallel to the simulated 
shock's normal direction, the fluid quantities only vary in 
the X direction, (ii) Collect the moments. Summation of 
particle masses and velocities are collected on a background 
computational grid based on PIC techniques. In this sub- 
step, the statistical average bulk speed of each grid repre- 
sents the velocity of each scattering center. Once the value 
of the bulk speed drops to zero, the position of the shock 
front is decided by the displacement of the corresponding 
grid, and it means that the shock position is moved with an 
evolutional velocity Vsh far away to the stationary reflected 
wall. Simultaneously, the size of the downstream region is 
extended dynamically with a constant velocity Vsh- Simi- 
larly, the forcshock region or precursor with a bulk velocity 
gradient is formed by the "back pressure" of the backward 
diffused particles. The moving of the FEB is also parallel to 
the shock moving with the same constant velocity Vsh- (hi) 
Applying multiple scattering laws. According to the scat- 
tering rate (i.e. Rs = dt/r, where Rs is the probability of 
the scattering events in time step dt, and r is the average 
scattering time). These fraction of the particles axe cho- 
sen to scatter the background scattering centers with their 
corresponding scattering angles obeying to the given Gaus- 
sian distributions. The chosen particles scatter off the col- 
lected background with their local velocities and scattering 
angles. The scattered particles move along their path until 
they have new scatters. In the duration of the time step, 
if the all chosen particles have completed their scatters, the 
background bulk speed is subsequently changed. In the turn, 
the varied background bulk speed also will change the par- 
ticle's individual velocity in the local frame in the next time 
step. The entire simulation time consists of the number of 
(Nt = 72000) time step involving the above three substeps. 

These presented simulations are all based on one- 
dimensional simulation box and the all simulated parame- 
ters has been described in detail elsewhere [Wang and Yan, 
2011). Here we list the simulation parameters in the Table 
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Figure 1. Schematic diagram of the simulation box. The shock is produced by incoming flow toward 
the reflective wall at the right boundary of the box (Xmax=300, Radial distance R=50). 



1. Upstream supersonic flow Uo with an initial Maxwellian 
thermal velocity Vl in their local frame and the inflow 
in a "pre-inflow box" (PIB) are both moving along one- 
dimensional simulation box from the left to the right. The 
parallel magnetic field Bo is along the x axis direction. FEB 
with a constant length Xfet = 90 in front of the shock 
position. The radial reflective boundary (RRB) is set as 
Ry,z = 50. The simulation box is dynamically consist of 
three regions: upstream, precursor and downstream. The 
bulk fluid speed in upstream region is U = Uo, the bulk 
fluid speed in downstream region is U — 0, and the bulk 
fluid speed with a gradient of velocity in the precursor re- 
gion is Uo > U > 0. To obtain the detailed information of 
the total particles in the simulation processes at any instant 
of time, we should build a large database for recording the 
velocities, positions, and the elapsed time of the all parti- 
cles, as well as the indices and the bulk speeds of the total 
grids. Then we can obtain the energy spectrums from the 
downstream, precursor, and upstream regions. The escaped 
particles' mass, momentum, and energy losses via the FEB 
can be also obtained. By analyzing the particle injection 
in the downstream region and the energy losses via FEB in 
the precursor region, we can find that how the prescribed 
scattering law to affect the shock compression ratio and the 
energy spectral index. 

To examine the relationships between the shock energy 
spectral index and the prescribed scattering law by the en- 
ergy injection and loss functions, we perform the Monte 
Carlo simulations with multiple scattering angular distri- 
butions using a new simulation system based on Matlab 
platform. The simulated cases are presented by Gaussian 
function with a standard deviation a and an average value 
(i.e., the expect value) jj, — involving four cases: (1) Case 
A: a = 7r/4. (2) Case B: ct = tv/2. (3) Case C: a = tt. (4) 
Case D: isotropic distribution. 



Table 1. The Simulation Parameters 



Physical parameters 


Dimensionless Value 


Scaled Value 


Inflow velocity 


Mo = 0.3 


403km/s 


Thermal speed 


DO = 0.02 


26.9km/s 


Scattering time 


T = 0.833 


0.13s 


Box size 


^max = 300 


lORe 


Total time 


tmax = 2400 


6.3minutes 


Time step size 


dt = 1/30 


0.0053s 


Number of zones 


nx = 600 




Initial particles per cell 


no = 650 




FEB distance 


Vf,b = 90 




Radial distance 


Ry,z = 50 


~ 1.5Re 



Note: The Mach number M =11.6. The Re is the Earth's radii. 
The data adapt from the Earth bow shock [Knerr et. al, 1996]. 



3. Energy analysis 
3.1. Shock structures 

We present the entire shock evolution with the velocity 
profiles of the time sequences in each case as shown in Fig- 
ure 2. The continuous infiow with a supersonic velocity Uo 
move from the left boundary {X = 0) of the upstream region 
to the downstream region at the right of the box with the 
time. The total bulk speed profiles are consist of three re- 
gions with the time: the upstream region U = Uo , precursor 
region < U < Uo, and downstream region U — 0. Total 
profiles of the bulk speed is distinct by two positions of the 
FEB and the shock front with the time. From the Cases A, 
B, and C to D, the precursor explicitly shows an increasing 
slope of the bulk speed, the shock's position Xsh also shows 
an increasing displacement increment in the X axis at the 
end of the simulation, respectively. This means the shock 
evolutes with an increasing velocity Vsh from the Cases A, 
B, and C to D, respectively. The simulated results of the 
dynamical shock in four cases are listed in the Table 2. In 
addition, by introducing a radial reflective boundary (RRB) 
in the present simulations, we also obtain the difference of 
the shock front position AXsh compared with the previ- 
ous simulations [Wang and Yan, 2011] with an increasing 
value of the {AXsh)A=-3.5, {AX,h)B=+6, {AX,h)c=+6, 
and {AXsh)D=+17.5 from the Cases A, B, and C to D, re- 
spectively. It is obvious to see that the affection of the RRB 
enhances this difference of the simulated shock for the four 
cases using the multiple scattering angular distributions. 

According to the relationships between the upstream and 
the downstream, we are able to calculate the total shock 
compression ratio rtot in the shock frame in each case as 
foUowings. 

Uo + \Vsh\ 
rtot = — pfT-j — (3) 

\Vsh\ 

Different shock evolutional velocity Vsh in different cases 
will probably lead to different dynamical shock structure. 
To showing this difference, we present the subtle velocity 
profiles at the end of simulation in each case in Figure 3. 
Evidently, the fluctuation of the velocity between the Vsub 
and Vd with an obliviously increasing value from the Cases 
A, B, and C to D, respectively. And the specific structure 
in each plot consists of three main parts: precursor, sub- 
shock and downstream. The smooth precursor with a large 
scale is between the FEB and the subshock's position Xaub, 
where the bulk velocity gradually drops from Uo to Vsub- 
The sharp subshock with a short scale just spans three-grid- 
length involving a deep drop of the bulk speed abruptly from 
Vaub to Vd, where the scale of the three-grid-length is about 
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Case A (o-i/4,|ii-0 ) 



Case B (o-i/2,|ii-0 ) 





Figure 2. The entire evolutional velocity profiles in four cases. The dashed line denotes the FEB posi- 
tion in each plot. The precursor is located in the area between the downstream region and the upstream 
region in each case. 
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Figure 3. The subtle structures of the subshock in four cases at the end of simulation time. The drops 
of velocity in the subshcok region are denoted by the values between Vsnb and Vd in each case. 



the thermal mean free path of the thermalized particles in 
the downstream region. So the subshock's velocity can be 
defined by the value of the Vsub in each case. The velocity 
Vd represents the downstream bulk speed at the shock po- 
sition at the end of the simulation. The bottom solid line 



denotes the backward shock evolutional velocity Vsh with 
an increasing value from the Cases A, B, and C to D, re- 
spectively. Because the subshock is the fraction of the total 
shock, we can calculate the subshock's compression ratio 
Tsub according to the total shock compression ratio rtot as 
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Figure 5. The four plots denote the mass losses, momentum losses, energy losses via the FEB and the 
injected energies in the downstream region, respectively. The solid line, dashed line, dash-dotted line and 
the dotted line represent the cases A, B, C and D in each plot, respectively. The units are normalized to 
the initial box proton mass Mp, initial box momentum Po and initial box energy Eq, respectively. 



following. 



Tsub = -rr- ^ ^tot (4) 
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3.2. Energy injection &i losses 

We have monitored the energy of the total particles over 
the time in different regions with respect to all cases. Fig- 
ure 4 shows all the types of energy functions with time. The 
Etot is the energy summation of the total particles in the to- 
tal simulation system over the time. The Etox is the energy 
summation of the actual particles in the simulation box over 
the time. The Epn is the energy summation of the contin- 
uous new particles enter into the simulation box from the 
"preinflow box" over the time. The Edowi is the energy sum- 
mation of the all particles in the downstream region over the 
time. The Edow2 is the energy summation of the all particles 
which their local velocity over the value of the initial velocity 
Uo in the downstream region over the time. The EdowS is the 
energy Edow2 mirms the Ei„j, which is the initial individual 
particle's energy (i.e. £k = l/2m(7(? + 1/2 mvf)) summation 
of the injected particles from the "thermal pool" at the lo- 
cal velocity of Vl = Uo to the superthermal particles in the 
downstream region, over the time. The Efeb is the energy 
summation of the total particles in the precursor region over 
the time. The Eout is the energy summation of the all parti- 
cles escaped from the FEB over the time. Clearly, the total 
energy Etot in the simulation system at any instant in time 
is not equal to the actual box energy Etox at any instant in 
time in each plot. It is evident from the real-time functions 
in Figure 4, the non-linear divergence between the curves for 
Ebox and Etot is produced with a decreasing value from the 
Cases A, B, and C to D, respectively. Also, the energy loss 
function Eout is produced with a decreasing value from the 
Cases A, B, and C to D, respectively. Simultaneously, the 
difference between the energy functions Edow2 and Edows 
shows an increasing energy injection Einj from the Cases A, 
B, and C to D, respectively. 

As shown in Table 3, all the listed results of the particle 
injection and losses in each case are calculated at the end 
of the simulation (i.e. Tmao^ =2400). The Mioss, Pioss and 
Eioss are the mass loss, the momentum loss and the energy 



Table 2. The results of the shock simulation 



Items 


Case A 


Case B 


Case C 


Case D 


^sh 


199.5 


165.5 


146 


106.5 


XpEB 


109.5 


75.5 


56 


16.5 




0.1075 


0.1460 


0.1757 


0.2525 


Vd 


-F0.0207 


-0.0024 


+0.0144 


+0.0045 


Vsh 


-0.0419 


-0.0560 


-0.0642 


-0.0806 


r-tot 


8.1642 


6.3532 


5.6753 


4.7209 


^sub 


2.9258 


3.0910 


3.3246 


3.9734 




0.7094 


0.7802 


0.8208 


0.9031 


^ sub 


1.2789 


1.2174 


1.1453 


1.0045 




11.4115 


14.2978 


17.2347 


21.6285 


ErrorBar 


+0.0017 


-0.0022 


+0.0014 


-0.0025 



Table 3. The results of the particle injection and losses 



Items 


Case A 


Case B 


Case C 


Case D 


Mloss 


1037 


338 


182 


38 


Ploss 


0.0352 


0.0189 


0.0123 


0.0025 


-^loss 


0.7468 


0.5861 


0.4397 


0.0904 


Etot 


3.3534 


3.4056 


3.3574 


3.4025 


Efeb 


0.8393 


0.5881 


0.5310 


0.3397 


^do-wl 


2.1451 


2.5612 


2.6359 


2.6903 


Einj 


0.1025 


0.1912 


0.2873 


0.3955 


Rinj 


4.78% 


7.47% 


10.90% 


14.70% 


^loss 


22.27% 


17.21% 


13.10% 


2.66% 



Notes; The units of mass, momentum, and energy are 
normalized to the proton mass Mp. initial total momentum Pq 
and initial box energy Eq, respectively. 



loss of the particles escaped via to the FEB, respectively. 
The Eff-h, Einj, Etot, and Edowi, with the unit of an initial 
box energy Eq, are all the energy values in their respective 
statistical volumes at the end of simulation. The R inj rep- 
resents the rate of the energy injection Einj with the total 
downstream energy Edowi at the end of simulation. And the 
Rloss represents the rate of the energy losses Eioss with the 
total energy in the system Etot at the end of the simulation. 
These correlations are presented as follows. 

Einj = Edow2 ~ EdowS (5) 

Rinj — Einj / Edowi (6) 
-Eioss = Eout (7) 

Rloss = Eout / Etot (8) 

For the comparison, the mass loss, momentum loss, en- 
ergy loss and energy injection functions with the time are 
calculated in Figure 5. Since the simulation system are 
based on the computational calculations, the existence of 
the energy losses is inevitable. Figure 5 show that the mass 
loss, momentum loss, and the energy loss functions with a 
decreasing value in any instant of time from the Cases A, 
B, and C to D, respectively. Among of theses loss func- 
tions, the energy loss function shows a decreasing values of 
(Eioss)a = 0.7468, {Eioss)b = 0.5861, (Eioss)c = 0.4397, 
and {Eioss)d ~ 0.0904 at the end of the sirrmlation from 
the Cases A, B, and C to D, respectively. On the contrary, 
the energy injection function show an increasing values of 
{E,nj)A = 0.1025, iEi„j)B = 0.1912, iE^n,)c = 0.2873, 
and {Einj)D = 0.3955 at the end of the simulation from 
the Cases A, B, and C to D, respectively. By of the ex- 
istence of the energy losses in the simulation system, the 
shock compression ratios are naturally affected according 
to the Rankinc-Hugoniot conditions. Therefore, the dif- 
ference of the energy losses or injection produced by the 
prescribed scattering angular distributions can directly af- 
fect all aspects of the simulated shock including the subtle 
shock structures, compression ratios, maximum energy par- 
ticles, and the energy spectrums, as well as other aspects. 
It is just this self-consistent injection mechanism and PIC 
techniques which allow the energy injection and loss func- 
tions to be obtained. So the further energy analysis for the 
diffusive shock acceleration could be done easily. 

3.3. McLximum energy 

We select some individual particles from the downstream 
region at the end of the simulation for obtaining the plots in 
the coordinates of the phase, space and time. The trajecto- 
ries of the selected particles arc shown in Figure 6. Among 
of these selected particles in each case, one of these tra- 
jectories clearly shows the fully acceleration processes of 
the maximum energy particle which undergoes the multi- 
ple crossings with the shock front. The maximum value 
of the local velocity marked in each plot shows an increas- 
ing values of (VL„,ax)A = 11.4115, {VL^ax)B = 14.2978, 
{VLmax)c = 17.2347, and {VLmax)D = 21.6285 from the 
Cases A, B, and C to D, respectively. And the correspond- 
ing statistical error of the local velocity in each case is 
listed in the Table 3. Conscqucnthr, The cutoff energy at 
the "power-law" tail in the energy spectrum is given with 
an increasing value of (iJmax)A = 1.23 McV, {Emax)B='^-^^ 
MeV, (£',no:r)c=2.80 MeV and {Emax)D='iAl MeV from 
the Cases A, B, and C to D, respectively. As for the escaped 
particles, owing to their energies are higher than the cutoff 
energy, they are not available in the system by of their es- 
caping via the FEB eventually. Since the FEB is a constant 
distance (i.e. Xfeb = 90) in front of the shock and maintains 
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Figure 6. The individual particles with their local velocities vs their positions with respect to time 
in each plot. The shaded area indicates the shock front, the solid line in the bottom plane denotes 
the position of the FEB in each case, respectively. Some irregular curves trace the individual particle's 
trajectories near the shock front with time. The maximum energy of accelerated particles in each case is 
marked with the value of the local velocity, respectively. 



the parallel moving of the shock front in each case, once an 
accelerated particle diffuse beyond the position of the FEB, 
this particle will be excluded from the system. The Table 3 
shows the numbers of the escaped particles at the end of the 
simulation with a decreasing mass losses of {Mioss)a = 1037, 
{Mioss)b = 338, {Mioss)c = 182, and {Mioss)d = 38 from 
the Cases A, B, and C to D, respectively. Also the en- 
ergy statistical data exhibit the energy loss rate with a de- 
creasing value of {Rioss)a = 22.27%, {Rioss)b = 17.21%, 
{Rioss)c = 13.10%, and {Rioss)d = 2.66% from the Cases 
A, B, and C to D, at the end of the simulation, respectively. 

Except for the maximum energy particles, there are also 
common energetic particles are shown in the plots with some 
of them obtained finite energy accelerations from the multi- 
ple crossings with the shock and some of them do not have 
additional energy gains owing to their lack of probability 
for crossing back into the precursor. If the cutoff energy of 
the simulation system is not effected by the prescribed scat- 
tering angular distribution, these maximum energy particles 
in different cases should be identical or at least be similar 
equal in the range of error bar. But the actual difference 
of the cutoff energy particles in different cases should be 
contributed by the different prescribed scattering angular 
distributions dominating the different energy injection. 

3.4. Heating, acceleration & spectrum 

As shown in Figure 7, the four energy spectrums with the 
"power-law" tails represent the four cases, averaged over the 
precursor region, at the end of simulation, respectively. The 
thin solid curve with a narrow peak is the initial Maxwellian 
distribution in the shock frame. The four extended energy 
spectrums are all consist of two very different parts: the low 
energy part and the high energy part. The low energy part 
in the left side of the initial spectrum, range from the low 
energy to the central peak, shows the "irregular fluctuation" 
in each case. The high energy part in the right side of the 



initial spectrum, range from the central peak energy to the 
cutoff energy, shows the smooth "power-law" tail in each 
case. The "irregular fluctuation" indicates that the super- 
sonic upstream fluid slows down in precursor region and its 
translational energy begin to convert into the irregular ran- 
dom energy. The "power-law" tail implies that the injected 
particles from the "thermal pool" in the downstream region 
scatter into the precursor region crossing the shock front 
for multiple energy gains and become into the superthermal 
particles. 

Look at extended curves closely, the low energy part in 
each case has a clearly joint point with the high energy part. 
And the joint point show an increasing energy value from 
the Cases A, B, and C to D, respectively. Consequently, 
the corresponding cutoff energy at the "power-law" tail in 
the precursor region also shows an increasing value from the 
Cases A, B, and C to D, respectively. This joint point should 
be correlated to the average thermal velocity in the down- 
stream region. As shown in the Figure 8, the four thermal 
velocity functions are averaged over the downstream region 
with the time. And each curve denotes the evolution of the 
average thermal velocity with the time and shows a constant 
after a certain duration (i.e, t = 500). Eventually, the av- 
erage thermal velocity Vth shows an increasing value from 
the Cases A, B, and C to D, at any instant of time, respec- 
tively. As expected, the energy injection from the "thermal 
pool" in the downstream region shows an increasing value 
from the Cases A, B, and C to D, respectively. Therefore, as 
show in the Figure 7, the energy spectrum in the precursor 
region shows an increasing hard spectral slope as the disper- 
sion value a of the Gaussian scattering angular distribution 
increases. This correlation of the energy spectrum averaged 
over the precursor region with the prescribed scattering law 
is consistent with the energy spectrum averaged over the 
downstream region. 
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Figure 7. This plot represents the energy spectrums on the precursor region at the end of the simulation. 
The thick solid line with a narrow peak at E =1.3105keV represents the initial Maxwell energy distri- 
butions. The solid, dashed, dash-dotted and dotted extended curves with the "power-law" tail present 
the energy spectrum corresponding to Cases A, B, C and D, respectively. All these energy spectrum are 
calculated in the same shock frame. 




Figure 8. This plot denotes the average thermal velocity with the time in the downstream region in each case. 
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Generally, we could predict the power-law energy 
spectral index from diffusive shock acceleration theory: 

dJ/dE oc E-^ (9) 

where dJ/dE is the energy flux and the F is the en- 
ergy spectral index. And the spectrum index can be 
calculated as following: 

Ttot = {rtot + 2)/[2x{rtot-l)]- (10) 

r«„6 = (r«„6 + 2)/[2x (r«„6-l)]. (11) 

According to Equation 10 and Equation 11, we substi- 
tute the corresponding values of the compression ratio 
r in each case. Then, the two types of energy spectral 
indices Ttot and Tgub in each case are calculated. As 
listed in the Table 2, the total shock energy spectral 
index shows an increasing value of the {rtot)A= 0.7094, 
iTtot)B =0.7802, (rtot)C =0.8208, and {T tot) D =0.9031 
from the Cases A, B, and C to D, respectively. How- 
ever, the subshock's energy spectral index is a decreas- 
ing value of the (r^„;,)^= 1.2789, {T^y_t,)B =1.2174, 
(rs„fc)C =1.1453, and (r,„fc)D=1.0045 from the Cases 

A, B, and C to D, respectively. 

As shown in Figure 9, all of the values of the sub- 
shock's energy spectral index are more than one (i.e. 
Tsub > 1 ); and the solid line denotes the subshock's 
energy spectral index with a decreasing value from the 
Cases A, B, and C to D as the energy injection in- 
creases, respectively. However, all of the values of the 
total shock's energy spectral index are less than one (i.e. 
^tot < 1), and the dashed line denotes the total shock's 
energy spectral index with an increasing value from the 
Cases A, B, and C to D as the energy injection increases, 
respectively. Simultaneously, as shown in Figure 10, the 
solid line denotes the subshock energy spectral index 
with a decreasing value from the Cases A, B, and C to 
D as the energy loss decreases, respectively. However, 
the dashed line denotes the total shock's energy spec- 
tral index with an increasing value from the Cases A, 

B, and C to D as the energy loss decreases, respectively. 
According to the diffusive shock acceleration theory, if 
the energy loss is limited to be the minimum, the simu- 
lation models based on the computer will more closely 
fit the realistic physical situation. The Figure 9 and 
Figure 10 indicate that the correlations of the energy 
spectral index with the energy injection or the energy 
losses are consistent with the energy spectral index is 
dependent on the prescribed multiple scattering angu- 
lar distributions. As seen from the Cases A, B, and C to 
D, the subshock's energy spectral index and the total 
shock's energy spectral index are both approximating 
to the realistic value one (i.e. F ~ 1 ) as the energy 
injection increases or as the energy loss decreases. As 
predicted, the Rankine-Hugoniot (RH) jump conditions 
allow to derive the relation of the compression ratio with 
the Mach number as: r = {-ja + l)/(7a - 1 + 



For a nonrelativistic shock, the adiabatic index 7^ = 
5/3 , if the Mach number M ^ 1, then the maximum 
compression ratio should be 4. According the Rankie- 
Hugoniot conditions, the total shock compression ra- 
tio should be less than standard value 4, and the cor- 
responding total shock's energy spectral index should 
be less than the standard value one for a nonrelativis- 
tic shock [Pelletier, 2001]. Simultaneously, we can see 
that if the energy injection achieves to the enough high 
level or the energy loss is limited to the enough low 
level, the subshock's energy spectral index will closely 
approximate the standard value of one. We present ex- 
plicitly these relationships between the energy spectral 
indices, the energy injection or energy losses, and the 
prescribed scattering law. And these relationships will 
be very helpful to improve simulation models by the 
best choice of the prescribed scattering law. Using the 
prescribed scattering law instead of the assumption of 
the transparent function in the thermal leakage mecha- 
nism, as far as the injection problem is concerned, the 
dynamical Monte Carlo model based on the PIC tech- 
niques is nothing less than a fully self-consistent and 
time-dependent model. 

4. Summary and conclusions 

In summary, we performed the dynamical Monte 
Carlo simulations using the Gaussian scattering angu- 
lar distributions based on the Matlab platform by mon- 
itoring the particle's mass, momentum and energy at 
any instant in time. The specific energy injection and 
loss functions with time are presented. We successfully 
examine the correlation between the energy spectral 
index and the prescribed Gaussian scattering angular 
distributions by the energy injection and loss functions 
in four cases. Simultaneously, this correlation is fur- 
ther enhanced by using the radial refiective boundary 
(RRB). 

In conclusion, the relationship between the energy 
injection or energy losses and the prescribed scatter- 
ing law verify that the shock energy spectral index is 
surely dependent on the prescribed scattering law. As 
expected, the maximum energy of accelerated particles 
is correlated with the particle injection rate from the 
"thermal pool" to superthermal population. So we find 
that the energy injection rate increases as the standard 
deviation value of the scattering angular distribution 
increases. In these multiple scattering angular distribu- 
tion scenario, the prescribed scattering law dominates 
the energy injection or the energy losses. So this self- 
consistent energy injection mechanism is capable to in- 
stead of the assumption of the thermal leakage injected 
function. Consequently, the cases applying anisotropic 
scattering angular distribution will produce a small en- 
ergy injection and large energy losses leading to a soft 
energy spectrum, the case applying isotropic scattering 
angular distribution will produce a large energy injec- 
tion and small energy losses leading to a hard energy 
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Spectral Index with Energy Injection 




Figure 9. The plot shows the correlation of the energy spectral index vs the energy injection. The 
triangles represent the total energy spectral index of the all cases. The circles indicate the subshock's 
energy spectral index of all cases. 



Spectral Index witti Energy Loss 
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Figure 10. The plot shows the correlation of the energy spectral index vs the energy losses. The 
triangles represent the total energy spectral index of the all cases. The circles indicate the subshock's 
energy spectral index of all cases. 



spectrum. These relationships will drive us to find a 
newly plausible prescribed scattering law which making 
the simulation model more close to the realistic physics. 
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